﻿using OSGeo.GDAL;
using OSGeo.OGR;
using OSGeo.OSR;
using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Threading.Tasks;

namespace ConsoleApp
{
    public class Gdal2tilesHelper1
    {
        public const string TILE_FORMAT = "png";
        public const string GRID_METADATA_JSON_FILE_NAME = "tilemap.json";
        private const int DEFAULT_PARALLELISM = 0;
        private static readonly object _reprojectLock = new object();

        public static void GenerateTiles(
            string inputPath,
            string outputDir,
            int minZoom = 0,
            int maxZoom = -1,
            int tileSize = 256,
            ResamplingMethod resampling = ResamplingMethod.Cubic,
            TileScheme scheme = TileScheme.XYZ,
            bool buildOverviews = true,
            int maxDegreeOfParallelism = DEFAULT_PARALLELISM,
            bool verbose = false)
        {
            if (!File.Exists(inputPath))
                throw new ArgumentException("输入文件不存在");

            if (Directory.Exists(outputDir))
                Directory.Delete(outputDir, true);
            Directory.CreateDirectory(outputDir);

            // ✅ 保留 GDAL 初始化
            GdalConfiguration.ConfigureGdal();
            GdalConfiguration.ConfigureOgr();
            Gdal.AllRegister();
            Ogr.RegisterAll();

            using (Dataset src = Gdal.Open(inputPath, Access.GA_ReadOnly))
            {
                if (src == null)
                    throw new InvalidOperationException($"无法读取文件: {inputPath}");

                string srcWkt = src.GetProjectionRef();
                SpatialReference srsSrc = new SpatialReference(null);
                if (!string.IsNullOrEmpty(srcWkt))
                    srsSrc.ImportFromWkt(ref srcWkt);

                SpatialReference srs3857 = new SpatialReference(null);
                srs3857.ImportFromEPSG(3857);
                string wkt3857;
                srs3857.ExportToWkt(out wkt3857, null);

                using (Dataset ds3857 = Gdal.AutoCreateWarpedVRT(src,
                                                                 string.IsNullOrWhiteSpace(srcWkt) ? null : srcWkt,
                                                                 wkt3857,
                                                                 GetGRA(resampling),
                                                                 0.0))
                {
                    if (ds3857 == null)
                        throw new InvalidOperationException("AutoCreateWarpedVRT 失败。");

                    if (maxZoom == -1)
                        maxZoom = CalculateDefaultMaxZoomForWebMercator(ds3857, tileSize);

                    if (verbose)
                        Console.WriteLine($"输入文件: {inputPath}, 输出目录: {outputDir}, 缩放级别: {minZoom}-{maxZoom}");

                    double[] gt = new double[6];
                    ds3857.GetGeoTransform(gt);
                    GetDatasetExtent(ds3857, gt, out double dsMinX, out double dsMinY, out double dsMaxX, out double dsMaxY);

                    var parOpt = new ParallelOptions
                    {
                        MaxDegreeOfParallelism = maxDegreeOfParallelism <= 0 ?
                                                 Environment.ProcessorCount : maxDegreeOfParallelism
                    };

                    for (int z = minZoom; z <= maxZoom; z++)
                    {
                        if (verbose) Console.WriteLine($"== 处理缩放级别 Z={z} ==");

                        TileRange tr = ComputeTileRangeForExtent(dsMinX, dsMinY, dsMaxX, dsMaxY, z, tileSize);
                        string zDir = Path.Combine(outputDir, z.ToString());
                        Directory.CreateDirectory(zDir);

                        Parallel.For(tr.MinX, tr.MaxX + 1, parOpt, tx =>
                        {
                            string xDir = Path.Combine(zDir, tx.ToString());
                            Directory.CreateDirectory(xDir);

                            for (int ty = tr.MinY; ty <= tr.MaxY; ty++)
                            {
                                int storeY = scheme == TileScheme.TMS ? ((1 << z) - 1 - ty) : ty;
                                string outPath = Path.Combine(xDir, $"{storeY}.{TILE_FORMAT}");
                                if (File.Exists(outPath)) continue;

                                GlobalMercator.TileBounds(tx, ty, z, tileSize,
                                    out double tMinX, out double tMinY, out double tMaxX, out double tMaxY);

                                try
                                {
                                    GenerateOneTile(ds3857, outPath, tMinX, tMinY, tMaxX, tMaxY, tileSize, resampling, verbose);
                                }
                                catch (Exception ex)
                                {
                                    if (verbose)
                                        Console.WriteLine($"瓦片生成失败 Z={z} X={tx} Y={ty} : {ex.Message}");
                                }
                            }
                        });
                    }

                    GenerateMetadataJson(ds3857, outputDir, minZoom, maxZoom, tileSize, scheme);
                }
            }
        }

        private static void GenerateOneTile(Dataset ds3857, string outFilePng,
            double tMinX, double tMinY, double tMaxX, double tMaxY,
            int tileSize, ResamplingMethod resampling, bool verbose)
        {
            OSGeo.GDAL.Driver mem = Gdal.GetDriverByName("MEM");
            OSGeo.GDAL.Driver png = Gdal.GetDriverByName("PNG");

            int srcBands = ds3857.RasterCount;
            int outBands = srcBands >= 3 ? Math.Min(4, srcBands) : 1;

            using (Dataset tileDs = mem.Create("", tileSize, tileSize, outBands, DataType.GDT_Byte, null))
            {
                // ✅ 设置投影和 GeoTransform
                tileDs.SetProjection(ds3857.GetProjectionRef());
                double[] gtTile = new double[] { tMinX, (tMaxX - tMinX) / tileSize, 0, tMaxY, 0, -(tMaxY - tMinY) / tileSize };
                tileDs.SetGeoTransform(gtTile);

                lock (_reprojectLock)
                {
                    CPLErr err = Gdal.ReprojectImage(ds3857, tileDs,
                        ds3857.GetProjectionRef(), ds3857.GetProjectionRef(),
                        GetGRA(resampling), 0.0, 0.0, null, null, null);
                    if (err != CPLErr.CE_None)
                        throw new InvalidOperationException("ReprojectImage 失败");
                }

                // ✅ 直接生成 PNG
                using (Dataset tmpCopy = png.CreateCopy(outFilePng, tileDs, 0, null, null, null))
                {
                    tmpCopy.FlushCache();
                }

                if (verbose)
                    Console.WriteLine($"生成瓦片: {Path.GetFileName(Path.GetDirectoryName(outFilePng))}/{Path.GetFileName(outFilePng)}");
            }
        }

        private static int CalculateDefaultMaxZoomForWebMercator(Dataset ds3857, int tileSize)
        {
            double[] gt = new double[6];
            ds3857.GetGeoTransform(gt);
            double res = Math.Max(Math.Abs(gt[1]), Math.Abs(gt[5]));
            int z = 0;
            for (int i = 0; i <= 30; i++)
                if (GlobalMercator.MetersPerPixel(i, tileSize) <= res)
                {
                    z = i; break;
                }
            return Math.Min(z + 1, 30);
        }

        private static void GetDatasetExtent(Dataset ds, double[] gt, out double minX, out double minY, out double maxX, out double maxY)
        {
            int w = ds.RasterXSize;
            int h = ds.RasterYSize;
            double[] xs = new double[4], ys = new double[4];
            PixToGeo(gt, 0, 0, out xs[0], out ys[0]);
            PixToGeo(gt, w, 0, out xs[1], out ys[1]);
            PixToGeo(gt, 0, h, out xs[2], out ys[2]);
            PixToGeo(gt, w, h, out xs[3], out ys[3]);
            minX = xs.Min(); maxX = xs.Max(); minY = ys.Min(); maxY = ys.Max();
        }

        private static void PixToGeo(double[] gt, double px, double py, out double gx, out double gy)
        {
            gx = gt[0] + px * gt[1] + py * gt[2];
            gy = gt[3] + px * gt[4] + py * gt[5];
        }

        private static TileRange ComputeTileRangeForExtent(double minX, double minY, double maxX, double maxY, int z, int tileSize)
        {
            GlobalMercator.MetersToTileXY(minX, maxY, z, tileSize, out int minTx, out int minTy);
            GlobalMercator.MetersToTileXY(maxX, minY, z, tileSize, out int maxTx, out int maxTy);
            int maxIndex = (1 << z) - 1;
            minTx = Math.Max(0, Math.Min(maxIndex, minTx));
            maxTx = Math.Max(0, Math.Min(maxIndex, maxTx));
            minTy = Math.Max(0, Math.Min(maxIndex, minTy));
            maxTy = Math.Max(0, Math.Min(maxIndex, maxTy));
            return new TileRange(minTx, minTy, maxTx, maxTy);
        }

        private static ResampleAlg GetGRA(ResamplingMethod m)
        {
            switch (m)
            {
                case ResamplingMethod.NearestNeighbor: return ResampleAlg.GRA_NearestNeighbour;
                case ResamplingMethod.Bilinear: return ResampleAlg.GRA_Bilinear;
                case ResamplingMethod.Cubic: return ResampleAlg.GRA_Cubic;
                case ResamplingMethod.CubicSpline: return ResampleAlg.GRA_CubicSpline;
                case ResamplingMethod.Lanczos: return ResampleAlg.GRA_Lanczos;
                default: return ResampleAlg.GRA_Cubic;
            }
        }

        private static void GenerateMetadataJson(Dataset ds3857, string outDir,
            int minZoom, int maxZoom, int tileSize, TileScheme scheme)
        {
            double[] gt = new double[6];
            ds3857.GetGeoTransform(gt);
            GetDatasetExtent(ds3857, gt, out double minX, out double minY, out double maxX, out double maxY);

            SpatialReference s3857 = new SpatialReference(null); s3857.ImportFromEPSG(3857);
            SpatialReference s4326 = new SpatialReference(null); s4326.ImportFromEPSG(4326);
            CoordinateTransformation ct = new CoordinateTransformation(s3857, s4326);

            double[] ll = new double[3] { minX, minY, 0 }; ct.TransformPoint(ll);
            double[] ur = new double[3] { maxX, maxY, 0 }; ct.TransformPoint(ur);
            double[] centerXY = new double[3] { (minX + maxX) * 0.5, (minY + maxY) * 0.5, 0 }; ct.TransformPoint(centerXY);

            var metadata = new
            {
                name = "gdal2tiles-csharp",
                description = "Tiles generated by C# Gdal2Tiles-compatible helper",
                version = "1.0.0",
                format = TILE_FORMAT,
                tileSize = tileSize,
                minZoom = minZoom,
                maxZoom = maxZoom,
                bounds = new[] { ll[0], ll[1], ur[0], ur[1] },
                center = new[] { centerXY[0], centerXY[1], (minZoom + maxZoom) / 2 }
            };

            string json = Newtonsoft.Json.JsonConvert.SerializeObject(metadata, Newtonsoft.Json.Formatting.Indented);
            File.WriteAllText(Path.Combine(outDir, GRID_METADATA_JSON_FILE_NAME), json);
        }
    }
    public struct TileRange { public int MinX, MinY, MaxX, MaxY; public TileRange(int minX, int minY, int maxX, int maxY) { MinX = minX; MinY = minY; MaxX = maxX; MaxY = maxY; } }

    public static class GlobalMercator
    {
        public const double EarthRadius = 6378137;
        public const double OriginShift = 2 * Math.PI * EarthRadius / 2.0;

        public static double MetersPerPixel(int zoom, int tileSize) => (2 * OriginShift) / (tileSize * Math.Pow(2, zoom));

        public static void MetersToTileXY(double mx, double my, int zoom, int tileSize, out int tx, out int ty)
        {
            double res = MetersPerPixel(zoom, tileSize);
            tx = (int)Math.Floor((mx + OriginShift) / (res * tileSize));
            ty = (int)Math.Floor((OriginShift - my) / (res * tileSize));
        }

        public static void TileBounds(int tx, int ty, int zoom, int tileSize,
                                      out double minX, out double minY, out double maxX, out double maxY)
        {
            double res = MetersPerPixel(zoom, tileSize);
            minX = tx * tileSize * res - OriginShift;
            maxX = (tx + 1) * tileSize * res - OriginShift;
            maxY = OriginShift - ty * tileSize * res;
            minY = OriginShift - (ty + 1) * tileSize * res;
        }

    }
}
